Inferencia bayesiana en la práctica: MCMC

datos <- read.csv(here::here("R_curso_modelos", "datos", 
                             "barbera_data_fire_total_climate.csv"))

# definimos parámetros de la previa
mu_a <- 0; sigma_a <- 1
mu_b <- 0; sigma_b <- 0.1

# Función que evalúa la densidad posterior no normalizada
# ("un" de "unnormalized", porque siempre colonializado, nunca incolonializado)
post_fire_un <- function(alpha, beta, log = T) { 
  
  ## Log Verosimilitud (mismo código que en like_fire)
  
  # Calculamos la media, lambda:
  lambda <- exp(alpha + beta * datos$fwi)
  # usando lambda evaluamos la verosimilitud de cada observación, 
  # en escala log:
  like_pointwise <- dpois(datos$fires, lambda, log = T)
  # la log-verosimilitud conjunta es la suma de las log-verosimilitudes por 
  # observación:
  like <- sum(like_pointwise)
  
  ## Log Previa
  lprior_a <- dnorm(alpha, mean = mu_a, sd = sigma_a, log = T)
  lprior_b <- dnorm(beta, mean = mu_b, sd = sigma_b, log = T)
  lprior <- lprior_a + lprior_b
  
  ## Log Posterior
  lpost <- like + lprior
  
  if (log) return(lpost)
  else return(exp(lpost))
}

# Creamos grilla de valores de alpha y beta.
side <- 200 # cantidad de valores por lado

alpha_seq <- seq(-0.75, 1.15, length.out = side)
beta_seq <- seq(0.1, 0.22, length.out = side)

grilla <- expand.grid(
  alpha = alpha_seq,
  beta = beta_seq
)

# Evaluamos la densidad posterior para cada valor
size <- side ^ 2 # cantidad de valores en la grilla
for (i in 1:size) {
  grilla$post[i] <- post_fire_un(grilla$alpha[i], grilla$beta[i], log = F)  
}

# Graficamos posterior no normalizada
ggplot(grilla, aes(x = alpha, y = beta, fill = post)) +
  geom_tile() + 
  scale_fill_viridis() +
  theme(legend.title = element_blank()) +
  xlab(expression(alpha)) + 
  ylab(expression(beta)) +
  ggtitle("Posterior no normalizada") +
  nice_theme()

size_a <- diff(alpha_seq)[1] 
size_b <- diff(beta_seq)[1]
area <- size_a * size_b
# área de cada celda en la grilla.

grilla$post <- normalize_dens(grilla$post, area)

# Graficamos posterior normalizada
ggplot(grilla, aes(x = alpha, y = beta, fill = post)) +
  geom_tile() + 
  scale_fill_viridis() +
  theme(legend.title = element_blank()) +
  xlab(expression(alpha)) + 
  ylab(expression(beta)) +
  ggtitle("Posterior normalizada") +
  nice_theme()

# Para facilitar las cuentas, pondremos los valores de la posterior en una
# matriz, análoga a la grilla (beta varía entre filas, alpha, entre columnas)
post_mat <- matrix(grilla$post, side, side, byrow = T)

alpha_tab <- data.frame(
  alpha = alpha_seq,
  post = NA
)

for (i in 1:side) {
  alpha_tab$post[i] <- sum(post_mat[, i] * size_b)
}
  
plot(post ~ alpha, alpha_tab, type = "l", xlab = expression(alpha),
     ylab = "Densidad posterior")

beta_tab <- data.frame(
  beta = beta_seq,
  post = NA
)

for (i in 1:side) {
  beta_tab$post[i] <- sum(post_mat[i, ] * size_a)
}
  
plot(post ~ beta, beta_tab, type = "l", xlab = expression(beta),
     ylab = "Densidad posterior")

joint <- 
ggplot(grilla, aes(x = alpha, y = beta, fill = post)) +
  geom_tile() + 
  scale_fill_viridis() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(legend.position = "none") +
  xlab(expression(alpha)) + 
  ylab(expression(beta)) +
  nice_theme()

ma <- ggplot(alpha_tab, aes(alpha, post, ymin = 0, ymax = post)) +
  geom_line(linewidth = 0.3) +
  geom_ribbon(color = NA, alpha = 0.1) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.line = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank())

mb <- ggplot(beta_tab, aes(beta, post, ymin = 0, ymax = post)) +
  geom_line(linewidth = 0.3) +
  geom_ribbon(color = NA, alpha = 0.1) +
  coord_flip() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.line = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank())

layout <- "
AA##
CCBB
"

ma + mb + joint + plot_layout(design = layout)

alpha_tab$mass <- alpha_tab$post * size_a # los pesos, suman 1
beta_tab$mass <- beta_tab$post * size_b

# Esperanza de alpha
ea <- sum(alpha_tab$alpha * alpha_tab$mass)

# Esperanza de beta
eb <- sum(beta_tab$beta * beta_tab$mass)

c(ea, eb)

# obtenemos la función de distribución acumulada marginal para alpha
alpha_tab$cdf <- cumsum(alpha_tab$mass)
plot(cdf ~ alpha, alpha_tab, type = "l", xlab = expression(alpha),
     ylab = "Probabilidad acumulada")

# creamos la función cuantil, es decir, la inversa de la acumulada.
# ajusta una spline. Al tomar como x la probabilidad acumulada (cdf) y como 
# Y los valores de alpha, es la inversa de la acumulada (que es cdf = f(alpha))
alpha_icdf <- splinefun(x = alpha_tab$cdf, y = alpha_tab$alpha, 
                        method = "monoH.FC")  
# esta función tiene como argumento la probabilidad acumulada (entre 0 y 1), 
# y devuelve el valor de alpha que acumula esa probabilidad (cuantil).
probs <- c(0, 0.025, 0.5, 0.975, 1)
quants <- alpha_icdf(probs)
names(quants) <- paste(probs * 100, "%")
quants

# (El mismo gráfico de antes)
p1 <- 
ggplot(grilla, aes(x = alpha, y = beta, fill = post)) +
  geom_tile() + 
  scale_fill_viridis() +
  theme(legend.title = element_blank()) +
  xlab(expression(alpha)) + 
  ylab(expression(beta)) +
  ggtitle("Posterior normalizada") +
  ylim(min(grilla$beta) - size_b, max(grilla$beta) + size_b) +
  xlim(min(grilla$alpha) - size_a, max(grilla$alpha) + size_a) +
  nice_theme()

# Remuestreamos celdas de la grilla con reemplazo, donde la probabilidad es 
# proporcional a la densidad posterior
nsim <- 10000 # número de muestras
rows_sim <- sample(1:size, size = nsim, replace = T, prob = grilla$post)

# Tomamos las filas que salieron sorteadas
draws_table <- grilla[rows_sim, ]

# Graficamos muestras
p2 <-
ggplot(draws_table, aes(x = alpha, y = beta)) + 
  geom_point(alpha = 0.1) + 
  ggtitle("Muestras de la posterior") +
  xlab(expression(alpha)) + 
  ylab(expression(beta)) +
  ylim(min(grilla$beta) - size_b, max(grilla$beta) + size_b) +
  xlim(min(grilla$alpha) - size_a, max(grilla$alpha) + size_a) +
  nice_theme()


(p1 | p2)

# Cada fila en draws_table es una muestra de la posterior conjunta (bivariada).
# Si queremos ver la marginal de alpha o beta, simplemente miramos esa columna
# en draws_table, ignorando los valores de la otra variable. La forma más rápida
# de resumir una marginal es graficando su densidad empírica.

par(mfrow = c(1, 2))
plot(density(draws_table$alpha), xlab = expression(alpha), main = NA)
plot(density(draws_table$beta), xlab = expression(beta), main = NA)

par(mfrow = c(1, 1))

# También podemos calcular medidas de resumen:
summary(draws_table[, c("alpha", "beta")])

# Para graficar ambas densidades en ggplot, alargamos:
draws_long <- pivot_longer(draws_table[, c("alpha", "beta")], # sólo esas columnas
                            cols = 1:2, names_to = "param", values_to = "value")

ggplot(draws_long, aes(x = value, color = param, fill = param)) +
  geom_density(alpha = 0.2) +
  facet_wrap(vars(param), scales = "free")

# Pr(alpha > 0)
sum(draws_table$alpha > 0) / nsim

# Pr(0 < alpha < 0.5)
sum(draws_table$alpha > 0 & draws_table$alpha < 0.5) / nsim

# Pr(exp(beta) > 1.2)
ebeta <- exp(draws_table$beta)
sum(ebeta > 1.2) / nsim 

# Pr(lambda > 20 | FWI = 18) // lambda es la media
FWI <- 18
lambda <- exp(draws_table$alpha + draws_table$beta * FWI)
sum(lambda > 20) / nsim

# Pr(y > 20 | FWI = 18) // y es la respuesta
FWI <- 18
lambda <- exp(draws_table$alpha + draws_table$beta * FWI)
ysim <- rpois(nsim, lambda)
sum(ysim > 20) / nsim

# Definimos la distribución de propuesta como una normal bivariada, y 
# aproximando la matriz de varianza-covarianza en base a la curvatura en el
# modo. Para ello, optimizamos y obtenemos la hessiana.

# optim necesita que la función a optimizar tome como argumento un vector 
# de parámetros
post_fire_opt <- function(x) {
  alpha <- x[1]
  beta <- x[2]
  lp <- post_fire_un(alpha, beta)
  return(lp)
}

opt <- optim(
  par = c(0 , 0),               # vector inicial 
  fn = post_fire_opt,           # función a optimizar
  method = "BFGS",              # método basado en gradientes (cool)
  control = list(fnscale = -1), # maximizar en vez de minimizar
  hessian = TRUE                # que calcule la hessiana
)

# La matriz de varianza-covarianza, asumiento normalidad es la inversa de la 
# -hessiana
V <- solve(-opt$hessian)

# Para muestrear con una normal multivariada de propuesta, V puede escalarse
# por un factor, que según la teoría, lo óptimo es 2.38 ^ 2 / d, con d siendo
# la dimensión de la posterior
optimal_factor <- 2.38 ^ 2 / 2

# Creamos una función que genera una cadena markoviana para nuestra posterior.
# start es un vector con el valor inicial de alpha y beta. n es el número de 
# muestras a tomar (incluye el start). factor es un factor que multiplica a 
# la matriz de varianza-covarianza de la propuesta; un buen valor es 2, pero es
# interesante variarlo para ver cómo afecta la eficiencia.
mcmc_fire <- function(start, n, factor = optimal_factor) {
  # Matriz para guardar las muestras
  draws <- matrix(NA, nrow = n, ncol = 2) 
  colnames(draws) <- c("alpha", "beta")
  
  # Vector para almacenar la log posterior, necesaria para comparar valores
  lp <- numeric(n)
  
  # Iniciamos la matriz de muestras y la lp
  draws[1, 1:2] <- start
  lp[1] <- post_fire_un(start[1], start[2], log = T)
  
  # Loop recursivo para hacer avanzar las coordenadas
  for (i in 2:n) {
    # Proponemos un nuevo valor (vector), simulando de la distribución de 
    # propuesta, centrada en draws[i-1, ] y con matriza de vcov V.
    prop <- rmvn(1, mu = draws[i-1, ], V * factor)
    
    # Evaluamos la probabilidad posterior en la propuesta
    lp_prop <- post_fire_un(prop[1], prop[2], log = T)
    
    # Calculamos el cociente de densidad posterior propuesta / valor anterior.
    # En escala log, la resta equivale al cociente:
    lp_diff <- lp_prop - lp[i-1]
    p_q <- exp(lp_diff) 
    # equivale a 
    # p_q = exp(lp_prop) / exp(lp[i-1])
    # pero en log es más estable
    
    # Aceptamos si p_q > 1 (mayor densidad en la propuesta), y sorteamos que 
    # la propuesta se acepte, con p_q si p_q es < 1.
    pkeep <- min(1, p_q)
    keep_prop <- ifelse(runif(1) < pkeep, T, F)
    
    if (keep_prop) { # si aceptamos la propuesta, la guardamos
      draws[i, ] <- prop
      lp[i] <- lp_prop   # para comparar con el siguiente
    } else {         # si no aceptamos, repetimos la muestra anterior
      draws[i, ] <- draws[i-1, ]
      lp[i] <- lp[i-1]   # para comparar con el siguiente
    }
  }
  
  return(draws)
}

# Para graficar sobre la grilla, la extendemos un poco. Esto permite ver cómo la 
# cadena se mueve por zonas de baja densidad también, probablemente porque la 
# iniciamos en una zona de baja densidad.
alpha_seq <- seq(-1, 1.5, length.out = side)  # límites más amplios que antes
beta_seq <- seq(0, 0.3, length.out = side)

size_a <- diff(alpha_seq)[1]
size_b <- diff(beta_seq)[1]

grilla <- expand.grid(
  alpha = alpha_seq,
  beta = beta_seq
)

for (i in 1:size) {
  grilla$post[i] <- post_fire_un(grilla$alpha[i], grilla$beta[i], log = F)
}

# Corremos una cadena corta comenzando dentro de la grilla
run1 <- mcmc_fire(start = c(-0.25, 0.15), n = 50)
draws_table <- as.data.frame(run1)

# Graficamos
ggplot(grilla, aes(x = alpha, y = beta, fill = post)) +
  geom_tile() + 
  scale_fill_viridis() +
  geom_point(data = draws_table, mapping = aes(x = alpha, y = beta),
             inherit.aes = F, size = 2, alpha = 0.3) +
  geom_path(data = draws_table, mapping = aes(x = alpha, y = beta),
             inherit.aes = F, alpha = 0.5, linewidth = 0.2) +
  theme(legend.position = "none") +
  xlab(expression(alpha)) + 
  ylab(expression(beta)) +
  ylim(min(grilla$beta) - size_b, max(grilla$beta) + size_b) +
  xlim(min(grilla$alpha) - size_a, max(grilla$alpha) + size_a) +
  nice_theme()

nsim <- 1000
set.seed(12359) # para reproducibilidad
run1 <- mcmc_fire(start = c(-0.5, 0.1), n = nsim)
draws_table <- as.data.frame(run1)

p1 <- 
ggplot(grilla, aes(x = alpha, y = beta, fill = post)) +
  geom_tile() + 
  scale_fill_viridis() +
  theme(legend.position = "none") +
  ggtitle("Densidad posterior") +
  xlab(expression(alpha)) + 
  ylab(expression(beta)) +
  ylim(min(grilla$beta) - size_b, max(grilla$beta) + size_b) +
  xlim(min(grilla$alpha) - size_a, max(grilla$alpha) + size_a) +
  nice_theme()

p2 <- 
ggplot(draws_table, aes(x = alpha, y = beta)) + 
  geom_path(linewidth = 0.1, alpha = 0.4) +
  geom_point(alpha = 0.1) + 
  ggtitle("Muestras") +
  xlab(expression(alpha)) + 
  ylab(expression(beta)) +
  ylim(min(grilla$beta) - size_b, max(grilla$beta) + size_b) +
  xlim(min(grilla$alpha) - size_a, max(grilla$alpha) + size_a) +
  nice_theme()

(p1 | p2)

draws_table$time <- 1:nsim

p3 <- 
  ggplot(draws_table, aes(time, alpha)) + 
  geom_line(linewidth = 0.3) +
  nice_theme() + 
  ylab(expression(alpha)) +
  xlab("Iteración")

p4 <- 
  ggplot(draws_table, aes(time, beta)) + 
  geom_line(linewidth = 0.3) +
  nice_theme() + 
  ylab(expression(beta)) +
  xlab("Iteración")

(p3 | p4)

nn <- 3000
run2 <- mcmc_fire(start = c(-2, 1), n = nn, factor = 0.1)
draws_table <- as.data.frame(run2)
draws_table$time <- 1:nn

p3 <- 
  ggplot(draws_table, aes(time, alpha)) + 
  geom_line(linewidth = 0.3) +
  nice_theme() + 
  ylab(expression(alpha)) +
  xlab("Iteración")

p4 <- 
  ggplot(draws_table, aes(time, beta)) + 
  geom_line(linewidth = 0.3) +
  nice_theme() + 
  ylab(expression(beta)) +
  xlab("Iteración")

(p3 | p4)

# Corremos 3 cadenas de igual largo comenzando en puntos al azar
iter_warmup <- 2000
iter_sampling <- 5000
n <- iter_warmup + iter_sampling
f <- optimal_factor * 0.05

c1 <- mcmc_fire(start = runif(2, -1, 1), n = n, factor = f)
c2 <- mcmc_fire(start = runif(2, -1, 1), n = n, factor = f)
c3 <- mcmc_fire(start = runif(2, -1, 1), n = n, factor = f)

# Las Combinamos en un array 3D, ya que son 3 matrices
arr <- abind::abind(list(c1, c2, c3), along = 3)
dimnames(arr) <- list(
  iteration = 1:n,
  variable = c("alpha", "beta"),
  chain = 1:3
)

# El paquete posterior requiere que las dimensiones sean iteration, chain, 
# variable, así que reordenamos:
arr <- aperm(arr, c(1, 3, 2))

# Converimos draws_array object, clase del paquete posterior
draws <- as_draws_array(arr)
summarize_draws(draws)

mcmc_trace(draws) # de bayesplot

# Quitamos el warmup, porque aún no había llegado al set típico (nos quedamos
# con las últimas iter_sampling de cada cadena)
iters_use <- (iter_warmup + 1) : n
draws_crop <- draws[iters_use, , ]
summarize_draws(draws_crop)

mcmc_trace(draws_crop)

functions {}           
data {}
transformed data {}
parameters {}              
transformed parameters {}
model {}                   
generated quantities {}

// (Esto se guarda en un archivo .stan)
data {
  int N;
  array[N] int y; // Número de incendios
  vector[N] x; // FWI
}

// Parámetros a muestrear
parameters {
  real alpha;
  real beta;
}

// Calculamos las cantidades derivadas
transformed parameters {
  vector[N] lambda = exp(alpha + beta * x);
}

// Acá definimos la log densidad posterior (o la que sea)
model {
  // Densidad previa
  alpha ~ normal(0, 1);
  beta ~ normal(0, 0.1);
  
  // Verosimilitud
  y ~ poisson(lambda);
}

model {
  // Densidad previa
  alpha ~ normal(0, 1);
  beta ~ normal(0, 0.1);
  
  // Verosimilitud
  y ~ poisson(lambda);
  
  /* 
    Esto equivale a decir "sumá a la log densidad posterior la 
    log densidad de alpha en una normal con media 0 y desvío 1, la
    log densidad de beta en una normal con media 0 y desvío 0.1, y la
    log verosimilitud (conjunta) de observar 'y' (es un vector) como realización
    de una Poisson cuyas medias son lambda (vector también).
  */ 
}

model{
  /*
    Stan guarda la log densidad objetivo en un objeto llamado 'target'.
    Cada vez que se evalúa la log posterior, target comienza valiendo 0. 
    Esta parte de la función se encarga de sumarle los términos correspondientes.
  */

  // Sumar log densidad previa
  target += normal_lpdf(alpha | 0, 1);
  target += normal_lpdf(beta | 0, 0.1);
  
  // Sumar log verosimilitud
  target += poisson_lpmf(y | lambda);
  
  // lpdf: Log Probability Density Function
  // lpmf: Log Probability Mass Function
}

model{
  /*
    Podemos mezclar ambas formas, pero siempre hay que usar una sola por 
    parámetro. (Si no, estará duplicada, lo cual es computacionalmente válido 
    pero sería sospechoso que tenga sentido.)
  */
  
  // Vale hacer
  alpha ~ normal(0, 1);
  target += normal_lpdf(beta | 0, 0.1);
  y ~ poisson(lambda);
  
  /*
    Las funciones de densidad o verosimilitudes suelen tener constantes
    normalizadoras. Usando target +=, Stan las conserva, pero usando los 
    sampling statements (~), los quita de la función, ahorrando un poco de 
    cómputo. En algunos casos necesitamos conservar esas constantes para 
    hacer una cuenta, y ahí conviene usar target +=.
  */
}

# El código de Stan, como escribimos arriba, debe estar guardado en un archivo
# .stan.

# Compilamos el modelo
model <- cmdstan_model(here::here("R_curso_modelos", "modelos", "nfuegos.stan"))

# Podemos revisar si lo escribimos bien
model$check_syntax()

# Creamos una lista nombrada para pasarle los datos. Los nombres de cada
# elemento de la lista tienen que ser exactamente los que definimos en 
# la sección data {}
stan_data <- list(
  N = nrow(datos),
  y = datos$fires,
  x = datos$fwi
)

# Y muestreamos la posterior
fit_mcmc <- model$sample(
  data = stan_data,
  chains = 4, 
  parallel_chains = 4,
  iter_warmup = 1000,
  iter_sampling = 1000
)

# Corremos los diagnósticos del muestreo. Si hay problemas, por defecto 
# aparecen warnings cuando sample() termina de correr, pero así podemos verlos 
# de nuevo.
fit_mcmc$cmdstan_diagnose()

# Resumimos las marginales usando el paquete 'posterior'
summ <- fit_mcmc$summary()
print(summ)

# Cuál es el peor neff y el peor rhat?
min(c(min(summ$ess_bulk), min(summ$ess_tail)))

max(summ$rhat)

# Seteamos un lindo color para bayesplot
color_scheme_set("viridisC")

# Extraemos las muestras y las formateamos como draws_array, una clase del
# paquete 'posterior'
draws <- fit_mcmc$draws(format = "draws_array")

# Gráfico de traza, para visualizar convergencia
mcmc_trace(draws, pars = c("alpha", "beta"))

# Gráfico de densidad por cadena
mcmc_dens_overlay(draws, pars = c("alpha", "beta"))

# Sólo para ilustrar, reajustamos el modelo pero guardando las iteraciones
# del warmup.
fit_mcmc_all <- model$sample(
  data = stan_data,
  chains = 4, 
  parallel_chains = 4,
  iter_warmup = 1000,
  iter_sampling = 1000,
  save_warmup = TRUE   # por defecto es FALSE
)

draws_all <- fit_mcmc_all$draws(format = "draws_array", inc_warmup = TRUE)
mcmc_trace(draws_all, pars = c("alpha", "beta"), n_warmup = 1000)

# Gráfico de dispersión, para ver un parámetro contra otro
mcmc_scatter(draws, pars = c("alpha", "beta"))

# Gráfico de pares, para ver de a muchos pares
mcmc_pairs(draws, pars = c("alpha", "beta"),
           off_diag_args = list(size = 0.75, alpha = 0.5))

# Avanzado: evaluar dóndo hubo transiciones divergentes
# Extraemos los parámetros del No-U-Turn Sampler (NUTS)
np <- nuts_params(fit_mcmc)
mcmc_pairs(draws, np = np, pars = c("alpha", "beta"),
           off_diag_args = list(size = 0.75, alpha = 0.5))

# (como no hubo transiciones divergentes, no muestra puntos rojos)

# Como vimos antes, podemos calcular probabilides como frecuencias.

# Extraemos las muestras como data.frame
d <- fit_mcmc$draws(variables = c("alpha", "beta"), 
                    format = "draws_df")
nsim <- nrow(d)

# Pr(exp(beta) > 1.2)
ebeta <- exp(d$beta)
sum(ebeta > 1.2) / nsim 

# Pr(lambda > 20 | FWI = 18) // lambda es la media!
FWI <- 18
lambda <- exp(d$alpha + d$beta * FWI)
sum(lambda > 20) / nsim

# Distribución posterior de lambda (la media) cuando el FWI es 18:
plot(density(lambda), xlab = expression(lambda), main = NA)

# Pr(y > 20 | FWI = 18) //
ysim <- rpois(nsim, lambda) # simulamos datos :)
sum(ysim > 20) / nsim

# Distribución predictiva posterior de y para FWI = 18:
plot(density(ysim, from = 0), xlab = "y", main = NA)

# Visualizamos ambos a la vez:
plot(density(lambda), xlab = expression(lambda, y), main = NA, xlim = c(0, 50))
lines(density(ysim, from = 0), xlab = "y", main = NA, col = 4)

# Y cuentas más arbitrarias aún. Ejemplo: ¿cuán probable es que el número de 
# incendios medio con un FWI = 20 sea el doble o mayor que cuando el FWI es 15?
# Esto requiere calcular la distribución posterior de lambda para FWI = 15 y 
# FWI = 20, luego calcular el cociente y evaluar qué proporción de muestras 
# tienen un cociente >= 2.
lambda_15 <- exp(d$alpha + d$beta * 15)
lambda_20 <- exp(d$alpha + d$beta * 20)
qlambdas <- lambda_20 / lambda_15
plot(density(qlambdas), main = NA, xlab = "Cociente de lambdas")
abline(v = 2, lty = 2)

# Pr(qlambdas >= 2)
sum(qlambdas >= 2) / nsim

# También podemos reponder lo mismo sobre el número observado de incendios crudo,
# no su media:
y_15 <- rpois(nsim, lambda_15)
y_20 <- rpois(nsim, lambda_20)
qy <- y_20 / y_15
plot(density(qy), main = NA, xlab = "Cociente de y")
abline(v = 2, lty = 2)

# Pr(qy >= 2)
sum(qy >= 2) / nsim

par(mfrow = c(1, 2))
plot(density(d$alpha), xlab = expression(alpha), main = NA, xlim = c(-3, 3))
curve(dnorm(x, 0, 1), add = T, col = 4)

plot(density(d$beta), xlab = expression(beta), main = NA, xlim = c(-0.5, 0.5))
curve(dnorm(x, 0, 0.1), add = T, col = 4)

par(mfrow = c(1, 1))